Synchronization by Nonlinear Frequency Pulling 
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We analyze a model for the synchronization of nonlinear oscillators due to reactive coupling and 
nonlinear frequency pulling motivated by the physics of arrays of nanoscale oscillators. We study the 
model for the mean field case of all-to-all coupling, deriving results for the onset of synchronization 
as the coupling or nonlinearity increase, and the fully locked state when all the oscillators evolve 
with the same frequency. 
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In the last decade we have witnessed exciting techno- 
logical advances in the fabrication of nanoelectromechan- 
ical systems (NEMS). Such systems are being developed 
for a host of nanotechnological applications, as well as for 
basic research in the mesoscopic physics of phonons and 
the general study of the behavior of mechanical degrees 
of freedom at the interface between the quantum and the 
classical worlds [l|, |2| • Among the outstanding features 
of nanomechanical resonating elements is the fact that at 
these dimensions their normal frequencies are extremely 
high — recently exceeding the IGHz mark [3j — facilitating 
the design of ultra-fast mechanical devices. Since with 
diminishing size output signals diminish as well, there 
is a need to use the coherent response in large arrays 
of coupled nanomechanical resonators (like the ones that 
were recently fabricated |j, |5J) for signal enhancement 
and noise reduction. One potential obstacle for achiev- 
ing such coherent response is the fundamental problem 
of the irreproducibility of NEMS devices. Clearly, as the 
size of a resonating beam or cantilever decreases to the 
point that its width is only that of a few dozen atoms, 
any misplaced atomic cluster dramatically can change the 
normal frequency or any other property of the resonator. 
Thus, it is almost inevitable that an array of nanome- 
chanical resonators will contain a distribution of normal 
frequencies. Here we propose to overcome this poten- 
tial difficulty by making use of another typical feature 
of nanomechanical resonators — their tendency to behave 
nonlinearly at even modest amplitudes. We shall demon- 
strate here that systems of coupled nonlinear nanome- 
chanical resonators (like the one we studied recently |6|) 
can self-synchronize to one common frequency through 



the dependence of their frequencies on the amplitude of 
oscillation. 

The synchronization of systems of coupled oscillators 
that have a distribution of individual frequencies is im- 
portant in many disciplines of science [7], |8J . The coher- 
ent oscillations can be used to enhance the sensitivity of 
detectors or the power output from sources, as proposed 
here. Synchronization is also important in biological phe- 
nomena, for example the collective behavior in popula- 
tions of animals, such as the synchronized flashing of fire 
flies, and the coherent oscillations observed in the brain. 

Although synchronization is often put forward as an 
example of the importance of understanding a nonlinear 
phenomenon, the intuition for the phenomenon, and in- 
deed the subsequent mathematical discussion, can often 
be developed in terms of simple linear ideas. Even the 
famous example of Huygens' clocks can largely be under- 
stood 9] in terms of a linear coupling of the two pendu- 
lums through the common support. The larger damping 
of the symmetric mode (coming from the larger, dissipa- 
tive motion of the support) compared with the antisym- 
metric mode tends to lead to a synchronized state with 
the two pendulums oscillating antiphase. The nonlinear- 
ity in the system is present in the individual motion of 
each pendulum in the mechanism to sustain the oscilla- 
tions, and to reach a full description of the range of pos- 
sible states must be included in the analysis. However, 
even without this drive the oscillators would still become 
synchronized through the faster decay of the even mode, 
albeit in a slowly decaying state. A second important fea- 
ture of the model describing the two pendulums, and of 
many other models used to show synchronization, is that 



the essential coupling between the oscillators is dissipa- 
tive, whereas in many physical situations the coupling is 
mainly reactive. 

In our example of the coupled array of nanomcchanical 
oscillators we expect to see predominantly reactive terms 
coming from the elastic coupling forces between the os- 
cillators. Furthermore, rather than the mode-dependent 
dissipation mechanism described above, we expect that 
for our nonlinear nanomechanical oscillators synchroniza- 
tion will arise from the intrinsically nonlinear effect of the 
frequency pulling of one oscillator by another. Thus, in 
this paper, we propose and analyze a model for synchro- 
nization involving reactive coupling between the oscilla- 
tors, which then leads to synchronization through non- 
linear frequency pulling — both effects must be present for 
synchronization to occur. 

Important advances in the understanding of synchro- 
nization have come from studying a simple model [lOj 
often known as the Kuramoto model |ll| . In this model, 
the oscillators are represented as phase variables, which 
in the absence of coupling simply advance at a rate that is 
constant in time, but with some dispersion of frequencies 
over the different oscillators. The coupling is included as 
infinite-range, or all-to-all coupling, so that the model is 
represented by the equations of motion for the N oscil- 
lators (with the dot denoting a time derivative) 
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Here the ojm are the individual oscillator frequencies 
taken from a distribution g{Lu), and K is a positive cou- 
pling constant. The synchronization is captured by a 
nonzero value of a complex order parameter ip 
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with the magnitude r,„ = 1 for the Kuramoto model. 

The Kuramoto equation shows rich behavior, includ- 
ing, in the large N limit, a sharp synchronization tran- 
sition at a value of the coupling constant K = K^ [ll|, 
which depends on the frequency distribution of the un- 
coupled oscillators g{uj). The transition is from an unsyn- 
chronized state with -(/i = in which the oscillators run at 
their individual frequencies, to a synchronized state with 
■0 7^ in which a finite fraction of the oscillators lock to 
a single frequency, Q — Q, equal to the mean frequency 
of the locked oscillators. The transition at Kc has many 
of the features of a second order phase transition, with 
universal power laws and critical slowing down |lll | , and 
a diverging response to an applied force [l3| . 

Equation |^ is an abstraction from the equations de- 
scribing most real oscillator systems, leaving out many 
important physical features. A natural generalization is 
to include the magnitude of the oscillations as dynamical 



variables [7| , while adding nonlinearities and considering 
reactive as well as dissipative coupling. Thus we are led 
to the model 
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The behavior including just nonlinear saturation and dis- 
sipative coupling (i.e. setting a = /3 = 0) was analyzed 
by Matthews et al. [l3|. We will instead study the case 
of reactive coupling (/? ^ 0, i^ = 0) and allow for non- 
linear frequency pulling (a ^ 0). This model then has 
two parameters: a the strength of the imaginary nonlin- 
ear term which yields the frequency pulling, and /? the 
reactive coupling strength. In addition the probability 
distribution g{u!) of the ujm must be specified. We will 
study the case of positive a and /3; for a symmetric dis- 
tribution g{u!) the results are the same changing the sign 
of both a and (3. 

The main focus of this paper is analyzing the behavior 
of Q, but first we want to show how such an equation 
might arise from the equations of motion of realistic non- 
linear coupled nanomechanical resonators. A possible set 
of equations describing such a system of N coupled res- 
onators (similar to the system we studied recently in a 
different context j.6|) is 
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D[Xm - k{Xm+l + Xm-l)] = 0. (4) 



The first two terms describe uncoupled harmonic oscil- 
lators, where the coordinate Xm measures the position 
of the m*'' nanomechanical cantilever or beam, oscillat- 
ing in its fundamental mode of vibration. We suppose 
the uncoupled oscillators have a linear frequency that is 
near unity (by an appropriate scaling of time) so that 
Sm <C 1. The third term is a negative linear damping, 
which represents some unspecified energy source to sus- 
tain the oscillations, and positive nonlinear damping, so 
that the oscillation amplitude saturates at a finite value. 
This saturation value is chosen to be of order unity by 
an appropriate scaling of the displacements Xm ■ The first 
three terms comprise a set of uncoupled van der Pohl os- 
cillators. The term ax^ is a reactive nonlinear term that 
leads to an amplitude dependent shift of the resonant fre- 
quency, observed experimentally in many nanomechani- 
cal resonators [ij, [ig • With u = this would give us 
a set of uncoupled Duffing oscillators. The final term 
is a nearest neighbor coupling between the oscillators, 
depending on the difference of the displacements. This 
is a reactive term, typical for either elastic or electro- 
static interaction between resonators that conserves the 
energy of the system. Others [lj| have considered non- 
linear oscillators coupled through their velocities; this is 
a dissipative coupling that would lead to X 7^ in the 
amplitude-phase reduction. 



The complex amplitude equation Q holds if the pa- 
rameters v^ a, D, 6m are sufficiently small. In this case 
the equations of motion are dominated by the terms de- 
scribing simple harmonic oscillators at frequency one. We 
may then write 
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where Zm (t) is slowly varying compared with the basic os- 
cillation frequency of unity, and • • • are correction terms. 
Substituting JHJ into the equations of motion Q and 
requiring that secular terms proportional to e'* vanish 
yields the amplitude equations 

2z,n ^ iiy + iSm)Zrn ^ [v + Szft) \Zm\ Zm 

- iD[Zrn - ^{Zm+l + Zm-l)]- (6) 

With a rescaling of time i — vt/2 © reduces to our 
model (jSJ except that in our model the nearest neighbor 
coupling is replaced by the all-to-all coupling convenient 
for theoretical analysis. 

Since we are interested in the behavior of Q for a 
large number of oscillators, it is convenient to go to a 
continuum description, where we label the oscillators by 
their uncoupled linear frequency uj = ojj rather than the 
index j, Zj -^ z[ijj). Introducing the order parameter 
(0) , the oscillator equations can be written in magnitude- 
phase form as 
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where ^ = 6* — is the oscillator phase relative to that 
of the order parameter, and lu is the bare oscillator fre- 
quency measured relative to A, which is the order param- 
eter frequency il — Q, shifted by — (a + /3) 
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At each time t the oscillators are specified by p{r, 9, ui, t), 
the distribution of oscillators at shifted frequency uJ over 
magnitude and phase values. The order parameter is 
given by the self-consistency condition 
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where g{Lo) is the distribution of oscillator frequencies 
expressed in terms of the shifted frequency (D. Note that 
unlike the cases of the Kuramoto model and Q with a = 
/3 = the imaginary part of this condition is not trivially 
satisfied even for the case of a symmetric distribution 
g{ui), and in fact serves to determine the frequency 11 of 
the order parameter. Furthermore, this frequency is not 
trivially related to the mean frequency of the oscillator 
distribution. 

To uncover more fully the behavior of our model Q we 
consider two issues: the existence of a fully locked state 



for large values of a/3; and the onset of synchronization, 
detected as the linear instability of the unsynchronized 
R = Q state. 

First we look at the fully locked solution for a bounded 
distribution of frequencies of width w. We define any 
state with an 0(1) magnitude of the order parameter R 
as synchronized. If all of the phases of a synchronized 
state are rotating at the order parameter frequency we 
call the state fully locked. The solutions are defined by 
setting dtr — Q which gives 
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and dfd = 0, which with Ijlll) can be written 

u, = F{9) = (3Rr^^{asm9-cos9), (12) 

where the solution to the cubic equation Hll() for r is to 
be used to form the function of phase alone F{9). The 
function F{9) acts as the force pinning the locked os- 
cillators to the order parameter. A particular oscillator, 
identified by its shifted frequency a), may be locked to the 
order parameter if (|12l) has a solution 9 = F~^{cu) and if 
this solution is stable. The stability is tested by lineariz- 
ing (|7I8(I about the solution. The fully locked solution 
will only exist if stable, locked solutions to (|12|l exist for 
all the oscillators in the distribution. In addition, the 
self consistency condition ifTUI) must be satisfied. 

For large values of a/3 the phases of the locked oscil- 
lators cover a narrow range of angles. The imaginary 
part of the self consistency condition p()|l shows that the 
range of phases must be around ^ = 0, and 112|l becomes 
(note r ~ 1 here) 
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The imaginary part of the self-consistency condition re- 
duces to (^) =0 (the average is over the distribution 
of frequencies), and the real part to simply _R ~ 1. Fi- 
nally, averaging H13|l over the distribution of frequencies 
fixes the order parameter frequency f2 ~ (lu) — a. This 
construction remains valid for large /3, so that unlike the 
case studied by Matthews et al. [13, "amphtude death" 
does not necessarily occur at large values of the coupling 
constant. The extension of this calculation to find the 
boundary of the fully locked state will be presented else- 
where. 

We now turn our attention to the initial onset of par- 
tial synchronization from the unsynchronized state. This 
is calculated by linearizing the distribution p around the 
unsynchronized distribution which is a uniform phase dis- 
tribution at r = 1, and seeking the parameter values at 
which deviations from the uniform phase distribution be- 
gin to grow exponentially. Care is needed in the analysis 
due to the important role the magnitude perturbations 
play. 

Introducing the small expansion parameter e charac- 
terizing the small deviations from the unsynchronized 



state, we write 

p{r,6,uj,t) ~ {2^r)-^6[r-l-eri{e,w,t)][l + eh{e,Ld,t)], 

as well as i? ~ £Ri, with ri,fi,Ri oc e'*'* with A the 
growth rate of the linear instability. With this expansion 
p remains normalized to linear order in e providing the 
average of /i over 9 is zero. The dynamical equations jT]- 
IH)) at 0(e) lead to the explicit solutions j'l = _Ri(Acos^+ 
B sin 9) with 
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The self-consistency condition p(J|) to first order in e gives 

l = iydLS.g(t:;)[(A + C)+*(i? + i?)], (14) 

We evaluate (I14|l at the onset of instability where the 
growth rate A — > (it is not sufficient to put A = since 
some terms of the integrals then diverge). We have eval- 
uated the integrals analytically for uniform, triangular, 
and Lorentzian distributions of bare frequencies. Here we 
present results for the triangular distribution, for which 
the resulting equation for the critical values of a, (3 and 
the order parameter frequency at onset must be solved 
numerically. 

Figure^ shows comprehensive results for a triangular 
distribution with full width w — 2. Panels (a-c) show 
the magnitude of the order parameter i? as a function 
of (3 for constant a scans from numerical simulations of 
(O for 1000 oscillators and K — 0. Limits of the unsyn- 
chronized state are consistent with the linear stability 
analysis. For the largest value shown, a = 0.9, the low 
(3 transition (3 = f3ci — 0.8 is weakly hysteretic, whereas 
the large (3 transition f3 < (3c2 = 3.7 is continuous. The 
state growing for (3 < /3c2 is a novel state with i? ^ 0, 
but with no oscillator frequency locked to the order pa- 
rameter, which has a frequency outside of the band of 
shifted oscillator frequencies. For f3 > 1.8 there is also a 
state with R close to unity in which all or most of the 
oscillators are locked to the order parameter. For smaller 



a = 0.42 there is a stable small R state for f3ci < f3 < f3c2, 
as well as a large R solution. For a = the large R syn- 
chronized state persists down to [3 > 1.6, whilst the un- 
synchronized state remains linearly stable for all /3 (panel 
c). Panel (d) shows the phase diagram, including results 
from the simulations as well as the linear stability anal- 
ysis of the unsynchronized and fully locked state. Over 
a large portion of the a, (3 plane there are two stable 
solutions — a large R synchronized state, and either the 
unsynchronized state (hatched region) or a small R state 
(cross hatched region) — leading to hysteresis for contin- 
uous parameter scans. Over the dotted portion only a 
synchronized state is stable, and over the unshaded re- 
gion only the unsynchronized state is stable. 
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FIG. 1: Results for a triangular distribution of full width w = 2. Panels (a)-(c) show the order parameter 7? found in numerical 
simulations of 1000 oscillators for sweeps increasing and decreasing /3 and for three representative values of a. The symbols 
are time-averaged values and the error bars are the standard deviations in R over the averaging time. Panel (d) shows the 
phase diagram deduced from sweeps at many values of q, and numerical calculations of the linear instabilities: solid line - 
linear instability of the unsynchronized state; short-dashed line - saddle-node line deduced from jumps of R in the numerical 
simulations (denoted by arrows in panels a-c); long dashed line - linear instability of the fully locked solution (the large R 
solution is fully locked to the right of this line). 



